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Abstract 

Spectral Energy Distribution (SED) fitting is a well-developed astrophysical tool 
that has recently been applied to high-redshift Lya-emitting galaxies. If rest-frame 
ultraviolet through near-infrared photometry is available, it allows the simultaneous 
determination of the star formation history and dust extinction of a galaxy. Lya- 
emitter SED fitting results from the literature find star formation rates ~ 3Mq yr^^, 
stellar masses ~ lO^M© for the general population but ~ IO^'^Mq for the subset 
detected by IRAC, and very low dust extinction, Ay < 0.3, although a couple of 
outlying analyses prefer significantly more dust and higher intrinsic star formation 
rates. A checklist of 14 critical choices that must be made when performing SED 
fitting is discussed. 
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1 Spectral Energy Distribution Fitting 



Although it is not widely recognized. Spectral Energy Distribution (SED) 
fitting and photometric redshift determination are identical. Both use the 
likelihood function L{z,T) that a galaxy's observed SED was generated by 
template type T at redshift z. This is typically generated from 
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i indexes the N photometry bands, p are template predictions in these bands, 
o are observed fluxes and a are observational errors. A best-fit normalization 
factor for each template/object combination can be solved for analytically to 
simplify optimization. 

When only a photometric redshift is desired, one marginalizes or optimizes 
over the nuisance parameter T, and it usually suffices to use a limited set 
of templates. SED fitting typically involves a large set of templates at fixed 
redshift, corresponding to a delta-function prior in z. This is justified when 
a spectroscopic redshift is available but is a dubious practice otherwise. A 
common practice is to determine photometric redshifts with a limited template 
set, fix the redshift at the best fit, and then perform SED fitting with a 
wide range of templates. This is computationally convenient but statistically 
inconsistent, is guaranteed to underestimate the uncertainties, and runs the 
risk of biasing the results significantly. 



Nonetheless, SED fitting is a robust, well-developed n iethod at low redshift, 
wher e spectroscopic redshifts are typically available (see lKannappan and Gawiserl . 



20071 and references therein for discussions of applications and caveats). Lya- 



Emitting (LAE) galaxies are well-suited to SED fitting because the narrow- 
band selection finds galaxies in a narrow enough redshift range that a fixed 
redshift can be assumed when fitting broad-band photometry. However, the 
signal-to-noise (S/N) available for photometry of these dim, high-redshift 
galaxies is considerably lower than for other galaxy types to which SED fitting 
has been applied. Hence we need to select SED methods that are appropri- 
ate for low S/N and to avoid over-fitting. Novel methods may be required to 
handle the unusual characteristics of these galaxies, in particular the guar- 
antee from high equivalent-width emission-line selection that a starburst is 
occurring at the time of observation, the corresponding guarantee that other 
nebular emission lines are strong, and the opportunity to utilize narrow-band 
photometry as part of the SED. 



Despite the low S/N available for LAEs, SED fitting has been claimed to 
allow the determination of star formation rate (SFR), stellar mass, stellar 
age, characteristic timescale for star formation (r), dust extinction, and a 
new q parameter descr ibing radiative transfer effects on Lya photons (e.g.. 



Finkelstein et al.l . 120071 ). When a starburst is occurring, the burst population 



is expected to dominate the rest-ultraviolet continuum. In the unusual case 
that dust extincti on is negligible, the SFR is proportional to the rest-UV flux 
( iKennicuttl . Il998l ) and the rest-UV slope depends on the age of the starburst 
and the shape of the initial mass function. In most galaxies, dust is too abun- 
dant for this interpretation to be meaningful, but in LAEs the selection met hod 
guarantees minimal dust extinction (see Fig. 14 of iGronwall et al.l . 120071 ). If 
dust extinction is not too high, the rest-NIR luminosity density is nearly di- 
rectly proportional to the stellar mass, with some dependence on the age of the 
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stellar population traced by rest- optical color since older main-sequence stars 
have a higher mass-to-light ratio (IBell et all l2003l : iPortinari et al.l . |2004J ). 



2 Results for Lya-emitting Galaxies 



Table [T] includes all reported results for emission-line-selected LAEs as of Oc- 
tober, 2008. This compilation does not include objects ini tially selected in 



the c ontinuum that turned out to show strong Lya emission (jPentericci et al. 



20091 ): these should be referred to as, e.g., Lya-emitting Lyman break galaxies, 
rather than LAEs. All observations find similar rest-UV fiux densities, hence 
analyses with large inferred dust extinction also report very large intrinsic 
SFR. Significant scatter exists in the age results, partially due to wide variation 
in model assumptions between constant SFR, r and instantaneous burst mod- 
els. The only robust determinations are of star formation rates ~ 3M0 yr~^, 
stellar masses ~ IO^Mq for the general population but ~ 10^'^Mq for the sub- 
set detected by IRAC, and of very low dust extinction. Ay < 0.3, although a 
couple of outlying analyses prefer significantly more dust and correspondingly 
higher intrinsic star formation rates. 

Discussion about the current results included concerns that systematic effects 
due to different analysis methods on these parameters can be serious. To 
estimate the severity of this, we agreed to trade SEDs to be analyzed by other 
groups using their methods to see how much difference is produced in the 
best-fit parameters and their reported uncertainties. Significant effort should 
be invested in developing a more consistent fitting method than the various 
ones used to obtain the current literature results. Determining evolution of 
LAE results with redshift or comparisons of LAEs and LBGs requires applying 
a common analysis method, as in the posters by S. Yuma and K. Ohta at this 
meeting. 



3 Discussion: 14 SED Fitting Choices 



In this section we discuss 14 choices necessary for any modeler to make when 
preparing to fit observed galaxy SEDs. 



1. Stellar population models: One typically choo ses a se t of stellar popu - 
lation models from BC03 fiBruzual and Charlotlliooih . M05 flMarastonl 120051 ) . 
or the unpublished CB08 (Charlot-Bruzual 2008) models. The M05 and CB08 
models include the empirically-determined contribution of thermally pulsat- 
ing asymptotic giant branch (TP-AGB) stars, which makes a major difference 
in rest-NIR photometry at intermediate ages. So BC03 should be avoided. 
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Table 1 
Results 



fro m F 09 (iFinkelstein et all \200^ . G 07 (ICawiser et all. [ioOTll. L07 



Lai et all lioOTl'l. L OS (|Lai et all \2004). NOT (jNilsson et all l2007l l. and P07 



([Pirzkal et al.l . |2007| ). NI (I) means that analysis was restricted to LAEs lacking 
(having) IRAC detections. Analyses by G07, LOS and NOT are of stacked popu- 
lations, whereas other results show the average and scatter amongst results for 
individual LAEs, with all error bars corresponding to the 6S% confidence level. 
Where two-population fitting was used (GOT, F09, POT), age and r refer to the 
younger population and young fraction is the fraction of the total stellar mass given 
in the 5th column that resides in the young population; if single-population fitting 
was used, young fraction is set to 1. If a constant SFR (instantaneous burst) was 
assumed, r is set to oo (0). The F09 results for radiative transfer of Lya photons 
yield q = l.Otll 
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although the difference is largest at solar metallicity and less severe for the 
l ower metallicities likely for LAEs. The discussion poin ted out that Starburst99 
( iLeitherer et al.l . Il999l : IVazquez and Leithererl . l2005l ) models include nebular 
emission, which the above models do not. They have been updated to in- 
clude TP-AGB stars (IVazquez and Leitheren . l2005l ) and should properly han- 
dle populations as young as 1 Myr, whereas the above models should not be 
trusted for young starbursts due to the lack of nebular continuum and emis- 
sion lines. However, the effects of metallicity and ionization parameter, both 
of which are unknown for LAEs, become important when using these models. 



2. Star formation history: In hierarchical cosmology, we expect galaxies 
to have complex star formation histories that combine a series of starbursts 
with quiescent star formation, with the average SFR increasing versus time 
during the early stages of galaxy formation before the supply of neutral gas 
becomes exhausted. In practice, SED-fitting approximates this using a smooth 
history, SFR(t) = SFR(to) exp{—(t — to)/T), where the SFR is assumed to 
be zero before to and the stellar age at the time of observation is tobs — to- 
A Constant Star Formation rate (CSF) corresponds to r ^ cx) and a Simple 
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Stellar Population (SSP) corresponds to a delta-function SFR at to (approx- 
imated by r — > 0). Negative values of r should be included in fits as they 
correspond to an exponentially increasing S FR which is what Av ould result 
from a constant specific star formation rate. iFinlator et al.l (120071 ) compared 
these smooth star formation histories with realistic bursty ones from cosmo- 
logical hydrodynamic simulations and found that both produced acceptable 
fits to high-redshift observations. 



3. Minimum age: The maximum age that a stellar population should be 
allowed to have is the age of the universe at the time of observation. This is 
mildly unrealistic, as a large SFR at the instant of the Big Bang makes little 
sense, but the observational consequences at z < 6 of beginning star formation 
at z = oo versus e.g., z = 15 are small given the short time difference. The 
minimum age is a subtler question. Following the discussion above, one should 
use Starburst99 (or the equivalent) for modeling any population younger than 
~ 10 Myr while varying the input metallicity and ionization parameter. 



4. Initial mass function (IMF): Although it is often stated that the 
IMF app ears universa l , a top-heavy IMF at high-redshift has recently been 
claimed (Ivan Dokkuml. l2008l : iDavel . |2008| ) and has been used to model LAEs 
( iLe Delliou et al.l . l2006l ) . Even different versions of the "universal" local IMF 
make a factor of two difference in stellar mass. This difference can be corrected 
for when comparing results as long as all papers report their assumed IMF. 
During the discussion it was suggested to use a Salpeter IMF between 100 and 
0.01 solar masses as a common standard to enable comparison of results. 



5. Metallicity range: Since the metallicities of LAEs have yet to be mea- 
sured, a conservative approach is to allow this parameter to vary between 
0.01 Zq and Zq. In theory, SED fitting can be used to fit the metallicity, 
but the effect is minor and is degenerate with other parameters, so it is best 
to consider this as a systematic uncertainty. The discussion pointed out that 
assuming a single metallicity causes an underestimate in the age uncertainty 
due to age-metallicity degeneracy in SED shape. 



6. Dust law: The Calzetti et al. (l2000f ) dust law gives the dust extinction 



as a function of wavelength, parametrized by E(B-V)= 0.3 Ay- Calibrated to 
local starburst galaxies, this is the dust law typically used for high-redshift 
SED fitting. However, we do not know if the dust law evolves with redshift, 
and given local variations we should consider this a significant uncertainty. It 
would therefore be useful to vary the dust law from Calzetti to SMC to Milky 
Way dust and to report the amplitude of uncertainties caused in other SED 
fitting parameters. The idea of removing dust from the fit was discussed; while 
defensible for LAEs where an absence of dust could be assumed, this could 
make it more difficult to compare results with other galaxy types where dust 
is clearly present. 



5 



7. IGM absorption: A model for absorption by tlie intergalactic medium 
(IGM) must be applied to template spectra. Most commo nly this f ollows 
Madaul (119951) but during t h e disc ussion updated analyses by ISongailal (120041 ) 
and iFaucher-Giguere et al.l (120081 ) were recommended. These references pro- 
vide a method for modeling IGM absorption of continuum emission, but the 
prescription that should be applied to Lya emission is unclear. Naively, all 
photons blue-wards of 1215. 67A rest-frame suffer IGM absorption, but the 
treatment in the literature varies from assuming no IGM effect on Lya to 
assuming that fully half the Lya photons have been absorbed by the IGM. 
Radiative transfer effects from resonant scattering and galac tic winds are too 



compl ex to yield an obvious recipe, as the simulations of IVerhamme et al. 



(12006! ) have shown. This provides a motivation for removing Lya emission 
from the SED fit by ignoring the narrow-band photometry, although more 
advanced options are described below. 



8. Number of stellar populations: Fitting a simple stellar population to 
an LAE is not a good approach if one is interested in determining age and 
total stellar mass. The presence of a starburst guarantees a young age and a 
low implied stellar mass for an SSP. Even CSF and r models are poorly suited 
for modeling a star formation history that is likely to have made a rapid 
upward jump at the beginning of the starburst (possibly due to a galaxy 
merger). Allowing exponentially increasing SFR through negative values of 
T will help somewh at. However, we nee d to utihze multiple - popul ation SED 
models, as done by iPirzkal et al.l (120071 ) and iGawiser et al.l (120071 ) . to reveal 
any underlying population of old stars whose star formation history does not 
smoothly tie onto the active starburst. 



However, as noted above, the low S/N available with LAE SEDs makes it 
important to avoid fitting too many parameters. Even an SED with 13 bands 
{U BVRIzJHK[3.6] [4.5] [5.8] [8.0]) can lead to degeneracies when fitting a dozen 
parameters plus the overall normalization. Each stellar population comes with 
4 degrees of freedom: instantaneous SFR, age, r, and Ay (stellar mass is de- 
termined by the first three). So we clearly cannot afford to fit more than 
two populations, and there are benefits to reducing the degrees of freedom 
by simplifying the populations e.g., by assurning a model at the age of the 
universe for the old population. iGawiser et al.l (120071 ) assumed that both pop- 
ulations see the same dust extinction; although this assumption is probably 
not true in general, it is very difficult to determine the dust reddening for an 
older population detected in only a couple IRAC bands where it dominates 
the photometry. It may make the most sense to combine a young model from 
Starburst99 with an old model from e.g. M05. 



9. Individual SEDs or stacked? Low S/N makes it difficult to obtain ro- 
bust results from SED fitting of individual LAEs, although several papers 
have chosen this approach. An alternative is to average ("stack") the pho- 
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tometry of an entire sample of LAEs. The latter approach achieves greater 
S/N at the cost of revealing "average" S EP parameters fo r the e ntire LAE 



population rather than individual objects. iFinkelstein et al.l (120091 ) compared 
these methods and demonstrated that the average of individual LAE SED 
parameters was somewhat similar to the SED parameters fit to the average 
LAE SED, although a more detailed investigation is needed to be able to 
properly compare stacked and individual SED fit parameters. The discussion 
pointed out that if an analysis splits a population of objects into subsets based 
upon, e.g., [3.6]/im flux, one should perform simulations to identify and sub- 
tract Malmquist-type biases caused by noise in the resulting SED parameters 
(noise causes the [3.6]/xm flux to be overestimated for the brighter stack and 
underestimated for the dimmer stack). 

10. Include nebular emission lines? An important decision in SED fit- 
ting is whether t o inclu de nebular emission lines in the template spectra. As 



Zackrisson et al.l (120081 ) have pointed out, these nebular emission lines can 
make a significant difference in SED fitting and photometric redshifts but are 
usually neglected. For very young populations at 4 < z < 5, Ha can make a 
significant contribution to the [3.6]yum photometry. BC03 models report the 
number of ionizing photons, which can be turned into Lya and Balmer series 
luminosities assuming Case B recombination. Starburst99 does this automat- 
ically and includes the critical [O II] and [O III] emission lines. 

11. Treatment of narrow-band photometry: Since spectral templates 
typically do not include emission lines and Lya is particularly complicated, 
the most common approach has been to ignore the narrow-band photometry 
of LAEs. Sometimes the inferred Lya flux is used to subtract the emission- 
line contribution from overlapping broad-band photometry. This is sensible, 
although the uncertainty in this correction should be propagated into the 
photometric uncertainties on the broadband fluxes. One caveat when using 
photometry from a narrow-band filter with a rounded (rather than top-hat) 
response curve is to avoid assuming that the LAEs are all at the redshift 
corresponding to the peak filter transmission, as this will systematically un- 
derestimate the objects' true Lya fiuxes. Ideally the redshift should be varied, 
and the filter transmission can at least be represented by a more typical value 



see 



Gronwall et al.l . 120071 ). 



A we l l-mot ivated attempt has recently been made by IFinkelstein et al.l (12007 



20081 . l2009l ) to incorporate Lya emission in the spectral templates to be com- 
pared with the full LAE SEDs including narrow-band photometry. These au- 
thors introduced a new SED parameter, g, where g = 1 implies trivial radiative 
transfer, g < 1 implies the expected preferential extinction of Lya photons, 
and q > 1 implies that Lya photons are enhanced versus continuum photons 
(ei ther by anisotr opic radiative transfer or possible clumpy dust as described 



by iNeufeldl . Il99ll ) . The combination of SFR, stellar age, and dust is sufficient 



7 



to predict the Lya flux when g = 1, and comparison with the continuum- 
subtracted narrow-band flux density determines q. The discussion appeared 
to produce agreement that in this sense q is independent of the broad-band 
SED fit and could be determined subsequently to simplify computation. 



12. Treatment of photometric uncertainties: When population-averaged 
fluxes are being used, a bootstrap analysis can be performed to include sample 
variance in the photometric uncertainties; this usually dominates the formal 
uncertainty in the average flux so is important to include. The discussion 
illuminated concerns about systematic errors in photometry, due primarily to 
the difficulty of performing aperture photometry in the IRAC bands given the 
much larger PSF and significant source confusion. 

Some statistically dubious habits have crept into the LAE SED fitting litera- 
ture with authors excluding "non-detections" from the fits or only penalizing 
the when a template flux exceeds a 3a upper limit. This alters the 
statistic to no longer follow a distribution, making interpretation difficult. 
Moreover, this modification is entirely unnecessary when SED fitting is done 
in flux (magnitudes should be avoided at all costs given the asymmetry of 
errors at low S/N). For low S/N data, noise fluctuations can cause negative 
fluxes, and these fluctuations are properly handled by feeding the observed 
fluxes and their formal photometric uncertainties. The related practice of only 
plotting photometry for bands with formal 3a detections is less dangerous but 
equally hard to justify. 



13. Method for determining best-fit model: Producing a likelihood for 
a single template is straightforward, but over a large parameter space the 
brute-force approach of determining the likelihood for every possible template 
can be very computationally intensive. Some authors are simplifying this using 
Markov-Chain Monte Carlo (MCMC), for which public routines exist. MCMC 
tries to avoid getting stuck in local minima, but extensive testing is still rec- 
ommended to be sure that the global optimum is being found successfully. 
Another decision t o be made is f requentist versus Bayesian interpretation of 



this likelihood (see iLuptonl . Il993l ). A frequentist analysis will simply call the 



template with maximum likelihood the best fit, whereas a Bayesian analy- 
sis will choose prior probabilities on all of the parameters and multiply the 
likelihood by these to report the best fit parameter set. 



14. Method for determining parameter uncertainties: The final choice 
to be made is how to determine the uncertainties in the best-fit SED pa- 
rameters. The most common method is to use Ax^ where the value is 
allowed to increase by the right amount to represent the 95% confidence level 
for the given number of degrees of freedom, ndof- However, frequentist anal- 
ysis actually recommends ruling out only those parameter values for which 
no template is a good fit to the data (despite considering all possible values 
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of other parameters) Tising absolute values of rather than Ax^. Unless the 
error bars are underestimated, this is quite conservative. Many astrophysicists 
instead use the Bayesian approach, where a "credible region" on a given pa- 
rameter is determined by projecting the likelihood function onto that single 
dimension, marginalizing over all other parameters, and keeping the parameter 
range that contains e.g. 95% of the posterior probability. To fully grasp pa- 
rameter dependencies and correlations, the full dimensional parameter space 
must be studied. Monte Carlo simulations can also be performed by varying 
the observed photometry within its reported uncertainties and determining 
the range of best-fit parameters. This is a reasonable approximation but is 
not statistically self-consistent. For example, a particular model might be an 
acceptable fit to all of these simulations but never end up as the best fit, so 
the parameter range of "best fits" can underestimate the true uncertainties. 
MCMC codes produce uncertainties along with their best fits. 

Discussion ensued over how to determine the number of degrees of freedom. 
It depends on whether a single model or an entire parameter space is being 
evaluated. For a model, ndof = ''T-data ~ '"-nuisance) where the final term is the 
number of nuisance parameters being fit such as the overall normalization. 
For a parameter space, ridof = JT-data — 'T-parameter, which essentially corresponds 
to treating all parameters as nuisance parameters to see if any model in this 
entire space is an acceptable fit to the data. 
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